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IBM 7094 PROGRAM FOR 
THE SIX-COIL PROBLEM 


INTRODUCTION 

One method of producing a uniform magnetic field with a very large volume 
of homogeneity is to use a six-circular-coil system. This report gives a brief 
description of the system, along with an analysis of its solution on the IBM 7094 
computer. 


SIX-COIL PROBLEM 


A theoretical discussion of the problem is given in the 1963 final report of 
the Goddard Summer Workshop Program. 1 

Our problem is to solve the following five equations: 


f 3 , i +y f 3 .2 + zf — =0 


3,3 


f 5,l + y S f 5,2 + zt f 5 , 3 =0 


f 7,l +y s2 f 7.2 + Zt2 f 7 , 3 =° Y 


( 1 ) 


f 9 , 1 +y s3 f 9,2 + Zt3 f 9 , 3 D ° 


: n,i + y g4 f 1 1 , 2 + zt4 f ii , 3 = 0 


J 

where the f's are the derivatives of Legendre Polynomials and are defined as 
follows: 


1 M. Speiser and D. L. Waidelich, "The Six-Circular-Coil System” in "Final Report of the Goddard 
Summer Workshop Program in Measurement and Simulation of the Space Environment” Publication 
X-320-63-264, Goddard Space Flight Center, Greenbelt, Maryland, pages C113-C121. 


1 



( 2 ) 


f, . = 5 x 2 - 1 

3.1 i 

f . = 21 x 4 - 14 x 2 + 1 

5.1 1 1 

f, . = 429 X 6 - 495 X 4 + 135 x 2 - 5 

7.1 i i i 

f . = 2431 x 8 - 4004 x 6 + 2002 x 4 - 308 x 2 + 7 

9.1 i i i i 

f, .= 29393 x 10 - 62985 x 8 + 46410 x 6 - 13650 x 4 + 1365 x 2 - 21, 

11.1 1 1 1 1 1 


where i = 1, 3 

This gives us five equations with seven unknowns: x i , x 2 , x 3 , y , z, s , and t. 
So that a finite number of solutions can be obtained, we shall assume values for 
x 2 and x 3> This gives us x 3 , y , z , s, and t to find. 

Initially we pick a trial value of x 1( and check to see if the first four equa- 
tions are satisfied. When two satisfactory trial values have been found, an at- 
tempt is made to satisfy the fifth equation by an iterative process that will be 
explained later. 

After a trial value of x t has been chosen (the method will be explained later) 
the first job is to find all the non-negative values of s that are no greater than 
some predetermined upper bound for s from the following formula. 


F(s) = a g6 s 6 + a 


65 SS +a 64 g4 +a 63 


s 3 + a g2 s 2 + a 61 S + a 60 = 0, 


61 


(3) 


where 


a 66 _ d 3 e 3 




a 65 = C 2 


(d 2 e 3 + d 3 e 2 ) 


64 


63 


- 2 Cj c 2 (dj e 3 + d 2 e 2 + d 3 e 3 ) 

(c 2 + 2 c 0 c 2 ) - (d Q e 3 + d 3 e 2 + d 2 e x + d 3 e Q ) 

2c o c i " < d o e 2 + d i e i + d 2 e o> 


> 


62 


a 61 = C 0 - ( d 0 e i + d l e 0> 
a 60 = " d 0 e 0 
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and where 


-f f f f f f -f f 2 f 2 f 

- r 3,l 5,2 r 5,3 1 7 , 1 7,3 *9,2 *3,3 *5,1 *7,2 l 9,. 




: 1 _ ^ *3,3 *5,1 *5,2 *7,1 *7,2 *9,3 *5,3 *7,3 ^ *3, 1 *5,2 *7,2 *9,1 

+ *3, 2 *5,1 *7,1 *9,2^ 


-f f f f f f - f f 2 f 2 f 

" 3,2 1 5 , 1 5,3 r 7,2 *7,3 9,1 3,3 *5,2 ‘7,1 


9,3 


> (5) 


^3 _ *7,2 *9,2 (*3,1 *5,3 *7,1 *3,3 *5,1 *7,3^ 

H-f f f f f f _f f 2 f 2 f 

2 “ r 3,3 1 5 , 1 5,2 7,1 7,3 9,2 3,1 ‘5,3 7,2 9,1 

^1 ~ *3,3 *5,1 *5,2 *7,2 *7,3 *9,1 " *3,2 *5,3 *7,1 *9,2 

^0 = *7,1 *9,1 ^ *3 , 2 *5,3 *7,2 ~ *3,3 *5,2 *7,3^ 


~\ 


(6) 


e 3 " *3,1 *5,1 (*5,2 *7,3 *9,2 *5,3 *7,2 *9,3^ 


~\ 


e 2 _ *3,1 

f 5,2 

3 

*7,1 

f 7,2 

*9.3 

f 3,2 

f 2 

5,1 

f 2 f 

1 7 , 3 1 9 , 2 

e i = *3,2 

*5,1 

f 5,3 

*7,1 

f 7,2 

f 9,3 

*3,1 

f 2 

5,2 

f 2 f 

1 7 , 3 *9,1 




e 0 " *3,2 *5,2 (*S,1 *7,3 *9,1 *5,3 *7,1 *7,3^ 


J 


(7) 


To solve our sixth-degree equation for s, we are using the Newton- Rapson 
Method, The basic formula for this method is 


S n + 1 = S " 


F (s ) 


" F'(s n ) 


( 8 ) 


3 


m 


where 


n 


F(s) 

F'(s) 


iteration number 

a 66 §6 + a 65 55 + a 64 34 + a 63 ^ + 3 62 S ' + 3 6 1 
6 a 66 S$ + 5 a 65 S “ + 4 V 53 + 3 a 63 ^ + Z 9 6 2 


S + a 60 
S + a 61 


Before the iteration is started, we must find an interval that contains one 
root. This is done by substituting values of s into F(s) until a change of sign 
occurs. We start with s = 0 and increment s by .01 until s equals some inter- 
mediate upper bound (predetermined) and then increment s by .1 until s equals 
the predetermined upper bound. If no change of sign occurs by the time s equals 
the upper bound our trial value of is discarded, and another trial value is used. 
If F(s)= 0 in one of the substitutions, the value of s is considered to be a root of 
F(s). In this program the intermediate upper bound has been set at s = 5, and 
the final upper bound at s = 10. If different values are desired, the change 
can be made by changing the source program and reassembling. 

Once two consecutive values of S that give opposite signs when substituted 
in F(s) are found, an initial estimate of the root must be made so that the itera- 
tion can be started. In making the estimate, we shall assume that F(s) is linear 
in the interval (a, b), so 


s _ F(a) (b - a) 
0 F(b) - F(a) 

Now, we go into the iterations where 


n + 1 


F ( s „> 

F '(0 


n = 0, 1, 2, 3, 


( 9 ) 


This iterative process is continued until 


As 

n 



< io - 8 
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If this is not accomplished in ten iterations, it is assumed that there is no solu- 
tion in this interval, and another interval is attempted. 

After a solution has been found or after it has been ascertained that there is 
no solution in the interval, an attempt is made to find another interval that con- 
tains a root. This process of seeking intervals and roots is continued until the 
upper bound on s is reached. 

When all of the s 's have been found, the next job is to find the values of t 
by the following formula 


t. 


B 

B 


1 , i 


2, i 


( 10 ) 




B I.i 


B 


2,1 


i.i 

f 3, 1 

+ f 7,l) 

A l,i 

f 3 , 3 

f 5 , 1 

2 , i 

f 5,l 

+ f 9,l) 

A 2 , i 

f 5,3 

f 7,l 


f 5.1 f 

f . f 


7,3 

A i,i 

f f 

3,3 1 5 , 1 

9,3 

A 2,l 

f f 

5,3 7,1 


> 




( 11 ) 


and where 


l.i 


5 i (^5,2 ^7,1 “ S X ^5,1 ^7,2 ) 
*3,2 ^5,1 - S i ^3,1 ^5,2 




> 


A, . = S i <f >.3 - 5 . i = 1, 2, • • • , < 6 

2,1 f f - s f f 

1 5 , 2 1 7 , 1 i 5,1 7,2 




( 12 ) 


Although all values of t are printed, only the positive values of t are used 
in the subsequent calculations. If there are no positive values of t , a note to 
this effect is printed, and the current trial value of is discarded as not being 
suitable (and another trail value of x x is substituted). 
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*A 


For each positive value of t that exists we shall find y and z from the follow- 
ing relations 



f - t f f 

3,3 r 5, 1 i 3, 1 A 5,3 


f - s f f 

3,2 r 5,l i 1 3 , 1 1 5 , 2 


(13) 


z 

1 


f 


3, 1 



r r 

1 3 , 2 + 1 3 , 3 


(14) 


y, 



i = 1, 2, 


< 6 


(15) 


We now have Xj , x 2 , x 3 , s , t , y, and z such that the first four equations are 
satisfied. Our next job is to satisfy the fifth equation-namely 


f 1 1 , 1 + y s4 f U ,2 + zt4 f 1 1 , 3 = 0 


To do this we shall use an iterative process of linear interpolation (or 
extrapolation). We compute the fifth equation as some function g(x x ) as shown 
below 


g( x !> = f n>1 + ys 4 f n>2 + zt 4 f n>3 


(16) 


We then seek another suitable x p and compute another g(Xj). The first estimate 
of the true x x is obtained by 


x ib g( x iA> - x ia g( x 1B > 
1 ” g(x 1A ) - g(x 1B ) 

where 


(17) 


x JA - first acceptance x x 
x ib second acceptance x 1 
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For the next approximation, x JB becomes x 1A , g(x JB ) becomes g( x 1A ), x x becomes 
x 1B , and another g(x t ) is computed to become the new g( x 1B ). This process is 
continued until 


1 - 


i,i 


1 , i + 1 


< 10 


_7 


i - iteration number 


or until a previously specified number of iterations has been made. When com- 
puting a new g(Xj) , new values of s, t, y, and z must be found because they are 
dependent uponxj . Also, since each x t can produce more than one s , and hence, 
more than one t , y, and z. There can be more than one g(x x ) for a given Xj ; 
therefore, more than one computed x x can be obtained from a pair of trial Xj 's. 
From this, we see that the process of keeping track of everything can become 
quite involved. 

When all of the final Xj's have been found, the last thing to do is to compute 
the values of s , t , y, and z for each x x . Also, in this last pass, four more sets 
of values are found for each x x ; these are 


r 2 = /s 


(18) 

r 3 =/I 


(19) 

u = y ( 

fl - x 2 \ 

| (20) 

2 r*> 

V 1 

z 

J q 

( l - 

| (21) 

3 

r 

1 3 

V-*lJ 


As a check these solutions are substituted in the left-hand side of the five 
equations, and the results are printed. 

At the start of this report I said that a couple values of x t are chosen, but I 
didn't say how they were obtained. There are several methods that can be used. 
In this program, lower bound is chosen, and this is the first trial value of . If 
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this value doesn't work, a given increment is added to this value. This process 
is continued until either a satisfactory x x is found, or until a predetermined 
upper limit is reached. With this x x added to the given increment being the 
lower bound, this process is repeated for the second x 1 . When the solutions 
have been obtained for the five equations, the entire process is repeated with 
the second acceptable trial being the lower bound. This goes on until the 
upper limit has been reached. 

This program was written for the IBM 7094 using Fortran IV. The appen- 
dices show the actual program along with the Input/Output specifications. With 
slight modification, this program should also work on other machines that use 
FORTRAN IV. 


INPUT 

The two input data cards are punched as follows: 


Card Number 

Columns 

Fortran symbol 

Format 

1 

1-14 

X2 

F14.7 


15-20 

blank 

- 


21-34 

X3 

F14.7 


35-40 

blank 

- 


41-42 

MAXRUN 

12 

2 

1-14 

XI 

F14.7 


15-20 

blank 

- 


21-34 

XMAX 

F14.7 


35-40 

blank 

— 


41-45 

DELTAX 

F5.3 


For a description of the FORTRAN symbols, see Appendix C. 


OPERATING AND MODIFYING INSTRUCTIONS 

This program is designed to be run under version 8 of the IBM 7094 IBSYS 
system using all seven index registers. 

As the program is written, the upper bound on the acceptable solutions of 
F(s) is s = 10. This value can be easily changed by changing the constant in the 
comparison in statement 320 (card number 11260) in the source deck and re- 
assemblying. 
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Appendix A 

FLOW CHART OF PROGRAM 
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FLOW CHART OF PROGRAM 












































Appendix B 
FORTRAN LISTING 
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no non 


$ IBFT C C0IL6 LIST. REF 10010 

DIMENSION X A ( 2 000 ) ,XB(2000) .XC( 2000) »GA<'2000) ,GB<2000) »XdUM(A) » 10020 

1F3«A) »F5< A) »F7< A) »F9( A) »F11 ( A) . S ( 7 ) . LT ( 2000 ) . A1 < 7 ) »A2 ( 7 ) »B1 ( 7 ) » 10030 

2B2<7) »T(7) »Y©Z(7) »Z( 7) .Y 1 7) »LTa< 2000) » R2 ( 7 ) » R3 < 7 ) » J2 ( 7 ) » J3 ( 7 ) » 100A0 

3LTB ( 2000 ) .CHECK (7.6) 10050 

WRUF (3.10) 10060 

10 F0RMATI 1H1///////////////////////////1H » AOX ♦ 16HSIX-C0I L PR0BLEM/11OO7O 

1H1 ) 10080 

IRUN=1 10090 

LMTA=1 10100 

LMTB=1 10110 

I SL1=0 10120 

I SL2=0 10130 

I SL3=0 101A0 

I SLA=0 10150 

I SL5=0 10160 

READ (2.20) X2.X3.MAXRUN 10170 

20 F0RMAT(2(F1A.7.6X),I2> 10180 

10190 

FIND ACCEPTABLE Xll AND X12 10200 

10210 

30 READ (2.A0) XI » XMAX , DELT AX 10220 

AO F0RMAT(2(F1A.7»6X) ,F5.3) 10230 

XB ( 1 ) =X1 10330 

WRITE(3.110)X2.X3 103A0 

110 F0RMAT(1H1////////////////1HO.3OX.61HC0MPUTaTI0N ©F SIX-C0IL PR0BL1O35O 

1EM WITH THE F0LL0WING INPUT DATA///1H .A0X.5HX2 = ,F1A.7//1|H .AOX, 10360 
25HX3 = , FI A . 7 ) 10370 

IPAG^=0 10380 

120 ICNTA=0 10A00 

ICNTB=0 10A10 

N A= 1 10A20 

NB= 1 10A30 

NC=1 10AA0 

130 I PAGE= I PAGE+1 10A50 

WRITE (3.1A0) IPAGE 10A60 

1A0 F0RMAT(1H1,11OX,5HPAGE ,IA) 10A70 

X I NT = XB ( NB > 10A75 

I F( ISL2-D1A7.1A7.1A3 10A76 

1 A3 WRITE (3.1A5) XINT 10A77 

1A5 F0RMATI 1H0.A3X.26HFINAL C0MPUTATI0N F©R Xl= .F1A.7///) 10A78 

G0 T© 155 10A79 

1A7 WRITE (3,150) XINT 10A80 

150 F©RMAT( 1H0.A0X.33HINTERMEDIATE C0MPUTAT I ®N F0R Xl= .F1A.7///) 10A90 

10500 

LEGENDRE P©LYN©MIALS 10510 
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C 10520 

155 WRITE (3*160) 10530 

160 F@RMAT(1H .20HLE6ENDRE POLYNOMIALS) 10540 

WRITE (3.170) 10550 

170 F0RMATI1HO.1HI »14X.5HF3( I ) .20X.5HF5I I ) »20X*5HF7 ( I ) »20X*5HF9( I ) .20X10560 
1 . 6f4Fl 1(1 )* 10570 

XDUM ( 1 )=XINT**2 10580 

XDUM ( 2 )=X2**2 10590 

XOUM ( 3 )=X3**2 10600 

D0 180 1=1.3 10610 

F3< I )=5.*XDUM( I)-l. 10620 

F 5 ( I ) = ( 2 1 • *XDUM ( I )-14. )*|XDUM( I ) + 1 . 10630 

F7( I ) = ( ( 429 • *XDUM ( I )-495. )*XDUM( I ) + 1 3 5 • )*XDUM( I ) — 5 • 10640 

F9( I )=( ( ( 2431 »*XDUM ( I >-4004. )*XDUM( I ) +2002. ) *XDUM( I )-308. )*XDUM( 1)10650 
1+7. 10651 

Fll ( I )=( ( ( (29393.*X0UM( I ) -6298 5. ) *XDUM ( I ) +46410. )*XDUM( I ) -13650. ) 10660 
1*XDUM ( I )+1365. )*XDUM( I ) — 2 1 . 10670 

180 WRITE (3.190) I . F3 ( I ) . F5 < I ) » F7 ( I ) » F9 ( I ) . F 1 1 ( I ) 10680 

190 F0RMAT(1H »I1»9X.4(E15.8.10X) .E15.8) 10690 

10700 

FIND C.D.E AND A F0R F(S) EQUATI0N 10710 

10720 

C2 = F3 < 1 >*F5<2 >*F5(3 )*F7( 1 ) *F7 ( 3 ) *F9 ( 2 ) -F3 ( 3 ) *F5 < 1 ) **2*F7 ( 2 > **2 107 30 

1* F9 ( 3 ) 10740 

0 = 2 ,*F3« 3 )*F5( 1 )*F5 (2)*F7( 1 )*F7( 2 )*F9 < 3)rF5 ( 3 > *F7 ( 3 ) * ( F3 < 1 ) *F5 ( 2 ) 1 0750 
1*F7(2 )j*jF9< 1)+F3(2)*F5< ll) *»F7 ( 1 )*F9(2) ) 10760 

C0 = F3 < 2 )*F5 ( 1 )*F5( 3)*F7(2 )*F7(3)*F9< 1 )-F3 ( 3 ) *F5 ( 2 ) ** 2 *F7 ( 1 ) **2 10770 

1* F9 ( 3 ) 10780 

D3 = F7 (2 ) *F9 ( 2 ) * ( F3 ( 1 ) *F5 < 3 ) **2*F7 ( 1 )-F3 <3 )*F5( 1 )**2*F7< 3) ) 10790 

D2 = F3 (3 )*F5 ( 1 )*F5( 2)*F7( 1 )*0 ( 3 ) *F9 ( 2 ) -F3 < 1 ) *F5 ( 3 ) **2*F7 ( 2 ) **2 10800 

1*F9 ( 1 ) 10810 

D1=F3 ( 3 )*F5 ( 1 )*F5<2 )*F7( 2)*F7(3)*F9( 1 )-F3 ( 2 ) *F5 ( 3 ) **2 *F7 ( 1 )**Z 10820 

1*F9(2) 10830 

D0=F7 ( 1 )*F9 ( 1 )*(F3(2)*F5(3)**2*F7(2)-F3(3)*F5<2)**2*F7<3) ) 10840 

E3=F3(1)*F5(1)*(F5(2)*F7(3)**2*F9(2)-F5(3)*F7(2)**2*F9(3) ) 10850 

E2=F3 ( 1 )*F5<2 )*F5(3)*F7( 1 )*F7(2)*F9(3)-F3<2)*F5< 1)**2*F7(3)**2 10860 

1*F9 ( 2 ) 10870 

E1=F3(2)*F5(1)*F5 ( 3)*F7( 1 ) *F7 ( 2 ) *F9 ( 3 ) -F3 < 1 >*F5 ( 2 ) **2*F7 ( 3 ) **2 10880 

1* F9 ( 1 ) 10890 

E0=F3 (2 )*F5(2 )*(F5( 1 ) *F7 ( 3 ) **2*F9 ( 1 ) -F5 ( 3 ) *F7 ( 1 )**2*F9 ( 3 ) ) 10900 

WRITE (3.200) D3»E3 10910 

200 F0RMAT ( 1H0.40X.5HD3 = »E 1 5 . 8 .2 OX , 5HE3 = .E15.8) 10920 

WR I TF (3.210) C2.D2.F2 10930 

210 F0RMATI1H » 5HC2 = .F15.8 .20X.5HD2 = ♦ F 1 5 . 8 » 20X . 5HE2 = »E15.8) 10940 

WR I TF (3.220) C1.D1.E1 10950 

220 F 0RMaT(1H ♦ 5HC 1 = . F 1 5 . 8 » 20X . 5HD1 = . FI 5 . 3 . 20X . 5HE1 = .E15.8) 10960 
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WRITE ( 3 » 230 ) CO » DO * E0 10970 

230 F0RMAH1H »5HC0 = , E 1 5 . 8 »2 OX , 5HD0 = » E 1 5 . 8 » 20X , 5HE0 = .E15.8) 10980 

A66=-D3*E3 10990 

A65=C2**2-(D2*E3+D3*E2) 11000 

A64 = 2.*Cl*C2-< D1*E3 + D2*E2+D3*E1 ) 11010 

A63=Cl**2+2»*C0*C2-( D0*E3+D1*E2+D2*E1+D3*E0 ) 11020 

A62=2.*C0*C1~( D0*E2+D1*E1+D2*E0 ) 11030 

A61=C0**2-(D0*E1+D1*E0) 11040 

A60=-D0*E0 11050 

WRITE (3*240) A66»A65»A64»A63»A62»A61»A60 11060 

240 F0RMAT ( 1 HO * 6X » 3HA66 » 1 6X , 3HA6 5 * 16X » 3HA64 * 16X , 3HA63 * 16X » 3HA62 » 16X * 11070 

13HA61 »16X,3HA60/lH * 6 ( E 1 5 . 8 * 4X ) , E 1 5 . 8 > 11080 

11090 

S0LUTI0NS 0F F(S) 11100 

11110 

A=0 11120 

1=1 11130 

FA=A60 11140 

IF(FA) 270,260*270 11150 

260 S( I ) =FA • 11160 

G0 T© 440 11170 

270 FD=FA 11180 

AO=A 11190 

A=A+ • 0 1 11200 

280 FA=( ( ( ( (A66*A+A65)*A+A64)*A+A63)*A+A62)*A+A61)*A+A60 11210 

IF(FA)300, 260*290 11220 

290 i F ( FD)380,310,310 11230 

300 IF(ED)310,310*380 11240 

310 IF(A-5. >270*320,320 11250 

320 IFIA-IO. >330,340*340 11260 

330 FD= FA 11270 

AD=A 11280 

A =A + • 1 11290 

G0 T© 280 11300 

340 I E ( 1-1 >350,350,370 11310 

350 IF(X1-XMAX)360, 353,353 11315 

353 WRITE <3*355) 11320 

355 F0RMAT ( 1H0,52HTHIS PRSBLEM CAN N®T BE S0LVED WITH THE GIVEN VALUES 1 1 32 5 
1///) 11330 

G® T© 1220 11335 

360 WRITE (3,365) XI 11340 

365 F0RMAT ( 1H0,4HX1= ,F14.7*28H IS N©T A SaTiSFACTSRY VALUE///) 11345 

X1=X1+DELTAX 11350 

XB ( 1 ) =X1 11353 

G© T© 130 11355 

370 LT(NB)=I-1 11360 
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G© T© 450 11375 

380 S( I )=-FD*(A-AD)/(FA-FD)+AD 11380 

I C= 1 11390 

390 H=( ( ( ( ( < A66*S< I )+A65 )*S( I )+A64 )*S( I ) +A63) *S( I ) +A62 )*S( I )+A61 ) *S ( I >11400 
1+A60 )/(((( ( 6 • *A66*S ( I ) +5 . * A65 ) *S ( I )+4.*A64)*S( I )+3»*A63)*S( \ ) +2 • 11410 

2*A62 ) *S ( I ) +A61 ) 11420 

S ( I ) =S ( I ) -H 11430 

I F ( ABS ( H ) ~ . 0 0000 00 1 ) 440 * 44 0 » 40 0 11440 

400 I F ( IC-10)410» 420*420 11450 

410 IC=IC+1 11460 

G© T© 390 11470 

420 WRlTF (3*430) S ( I ) . X A ( NA ) » XB < NB ) 11480 

430 F0RMAT(1HO*3HS= *F14.7,32H C0NVERGF5 T0© SL0WLY WHEN XH= .F14.7* 11490 
1 1 OH AND X12 = *F14. 7// ) 11500 

G® T© 310 11510 

440 1=1+1 11520 

G© T© 310 11530 

450 WRITE (3*460) 11540 

460 F0RMATI1HO.1HI *11X*4HS( I ) * 1 5 X , 5 HA1 ( I ) * 17X , 5HA2 ( I ) *17X»5HB1 ( I ) * 1 7X » 1 1 5 50 
15HB2( I )*17X*4HT( I) ) 11560 

LTD=LT(NB) 11570 

D© 470 1=1 *LTD 11580 

A 1 ( I ) =S ( I )*( F5<2)*F7( 1 )-S( I )*F5 ( 1 )*F7< 2 ) ) / < F3 < 2 ) *F5 ( 1 ) -S ( I )*F3< 1 ) 11590 
1 *F5 ( 2 ) ) 11600 

A 2 ( I ) =S( I )*( F7< 2 )*F9 < 1 )-S( I ) *F7 ( 1 )*F9 ( 2 ) ) / < F5 <2 )*F7 ( 1 ) -S ( I ) *F5 < 1 ) 11610 
1*F7 ( 2 ) ) 11620 

Bill) =F5 ( 3 ) * < A 1 ( I )*F3 ( 1 )+F7 ( 1 ) )*A2 < I ) *F5 ( 3 ) *F7 ( 1 )-F7 ( 3 )*( A2 < I ) 11630 

1*F5 ( 1 ) +F9 ( 1 ) ) *A1 ( I )*F3 ( 3 )*F5 (1 ) 11640 

B2 ( I ) =F5 ( 1 ) *F7 ( 3 )*A2 ( I )*F5 ( 3 )*F7( 1 )-F7( 1 ) *F9 (3)*A1( I )*F3(3)*F5(1) 11650 
T ( I ) =B1 ( I ) /B2 ( I ) 11660 

470 WRlTF (3*480) I * S ( I ) * A 1 < I ) * A2 ( I ) . B1 < I ) » B2 ( I ) » T < I ) 11670 

480 F©RMAT(1H *I1*7X*F11.7,5(7X*E15.8) ) 11680 

I F ( ISL2-1 ) 510,510,490 11690 

490 WRlTF (3*500) 11700 

500 F®RMAT (1H0*6X,4HT ( J) , 15X*4HY(J) »15X*4HZ(J) *14X*5HR2(J) . 14X , 5HR3 ( J ) 1 1 71 0 
1*14X,5HJ2( J) *14X,5HJ3( J) ) 11720 

G® T© 530 11730 

510 WRITE (3,520) 11740 

520 F©RMAT I1H0*1HJ* 14X»4HT( J ) *20X»6HY/Z( J ) *20X *4HZ( J ) *2lX *4HY (Z) *21 X, 11750 
1 4HG ( J ) ) 11760 

530 J=1 11770 

1=1 1 11780 

540 I F ( T ( I ) >550,550*610 11790 

550 I F ( I -LTD ) 560,570*570 11800 

560 1=1+1 11810 

G© T© 540 11820 
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570 IF( J-l )590»590.580 11830 

580 J=J-1 11840 

I F ( ISL2-1 >650,750,680 11850 

590 WRITE (3*600) 11860 

600 F0RMAT ( 1H0 » 23HN© POSITIVE T AVAILABLE) 11870 

I F ( ISL1-1 >350*350*603 11875 

603 IF(LMTB-1 1604*604,605 11880 

604 WR I TE ( 3 * 365 ) XINT 11883 

6© T© 713 11885 

605 IF(NB-LMTB)608*606*606 11887 

606 LMTB=LMTB-1 11890 

6© T© 770 11893 

608 LMTB=LMTB-1 11895 

D© 609 I =NB »LMTB 11897 

609 XB ( I ) = XB ( 1 + 1 ) 11900 

6© T© 130 11905 

610 Y©Z( J)=-(F3<3)*E5( 1 ) — T ( I )*F3 ( 1 )*F5<3) ) / ( F3 < 2 ) *F5 < 1 ) -S ( I ) *F3 ( 1 ) 11910 

1*E5 ( 2 ) ) 11920 

Z ( J ) =-F3 ( 1 )/<Y©Z(j)*F3<2)+F3(3) ) 11930 

Y( J)=Y©Z(J)*Z( J) 11940 

I F ( ISL2-1 1620,720*660 11950 

620 GA( J) =F11 ( 1 >+Y < J)*S( I )**4*F11 (2 ) +Z ( J )*T ( I )**4»F11 (3 ) 11960 

WRUF <3*630) J*T( i ) ,Y©Z(J) ,Z(J) ,Y(J) ,GA( J) 11970 

630 F©RMAT(1H *I1*9X,4(E15.8*10X).E15.8) 11980 

I F ( I-LTD)640, 650*650 11990 

640 J=J+1 12000 

1*1+1 12010 

G© T© 540 12020 

650 LTA ( NA ) = J 12030 

I SL2 = 1 12040 

I SL 1 = 1 12045 

653 IFIX1-XMAX)655*353*353 12050 

655 XA ( 1 ) =X1 12055 

X1=X1+DELTAX 12060 

XB ( 1 ) =X1 12063 

G© T© 130 12065 

660 R2( J)=SORT(S( I ) ) 12070 

R 3 ( J ) =SQRT ( T ( I ) ) 12080 

J 2 ( J)=Y( J)*(1.-X1**2 )/(R2( J)**3*(l.-X2**2) ) 12090 

J3( J) =Z( J)*( l.-Xl**2 ) /(R3( J) **3*( l.-X3**3) ) 12100 

WRITE (3*670) T ( I ) * Y( J ) ,Z ( J) ,R2 ( J) *R3 ( J ) * J2 < J ) * J3 ( J ) 12110 

670 E©RMAT(1H » 6 ( El 5. 8 * 4X ) , E 1 5 . 8 ) 12111 

CHECK ( J » 1 ) = F3< 1 )+Y( J)*F3 (2 )+Z( J)*F3( 3) 12112 

SDUM=S ( I ) *Y ( J ) 12114 

TDUM=T ( I ) *Z ( J ) 12116 

CHECK! J*2 )=F5( 1 )+SDUM*F5 (2)+TDUM*F5( 3 ) 12118 
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SDUM=SDUM*S( I ) 12120 

TDUM=TDUM*T( I ) 12122 

CHECK (J»3)=F7(1)+SDUM*F7(2)+TDUM*F7(3 ) 12124 

SDUM=SDUM*S( I ) 12126 

TDUM=TDUM*T< I) 12128 

CHECK! J*4)=F9 ( 1 )+SDUM*F9(2 ) +TDUM*F9 1 3 ) 12130 

SDUM=SDUM*S( I ) 12132 

TDUM=TDUM*T( I ) 12134 

CHECK! J»5 >=F11 < 1 )+SDUM*Fll ( 2 ) +TdUM*F 1 1 ( 3 ) 12136 

IF! I-LTD1640, 680,680 12138 

680 I F ( NB-LMT B)690»700»70Q 12140 

690 NB=NB+1 12150 

G© T© 130 12160 

700 WRlTF(3*703) <I*I=1*5) 12163 

703 F©RMAT!1H0»21HVALIDITY ©F S®LUTI©NS//1H , 1HJ , 5 < 8X , 9HEOUAT I ©N *11* 12165 

12 X ) ) 12767 

D© 705 I=1»J 12770 

705 WRlTF(3*707) I»<CHECK<I»IT)»IT=1,5) 12773 

707 F©RMAT(1H , II *51 5X*E15.8 ) ) 12175 

WRITE (3*710) X I NT »X2 » X3 12177 

710 F0RMATI 1H0.5HX1 = *F14.7/1H *5HX2 = »F14,7/1H *5.HX3 = *F14.7) 12180 

713 IF(X1-XMAX)715*1240*1240 12181 

715 I RUN = 1 12182 

LMTA=1 12184 

LMTB=1 12186 

I SL1 =0 12188 

I SL2=0 12190 

I SL3=0 12191 

ISL4=0 12192 

I SL5 = 0 12193 

XB ( 1 ) =X1 12194 

G© T® 120 12196 

720 JB=J+ICNTB 12200 

GB ( JB ) =F1 1 ( 1 ) +Y ( J ) *S ( I )**4*F11 < 2 ) +Z ( J ) *T ( I )**4*F11 (3 ) 12210 

WRITE (3*630) J * T { I ) , Y®Z ( J ) » Z ( J ) , Y ( J ) , GB ( JB ) 12220 

IF! I-LTD)730*750*750 12230 

730 IF(JB-2000)640. 740*740 12240 

740 WRITE (3*745) 12250 

745 F0RMAT ( 1H0.30HHELP - I AM BEING SQUEEZED ©UT///) 12260 

G® T© 1240 12270 

750 IF! ISL1-1 >753*753*755 12273 

753 I'SLl = 2 12275 

XB ( 1 ) =X1 12277 

755 LTB ( NB ) = J 12280 

IF(NB-LMTB)760, 770*770 12290 

760 NB=NB+1 12300 
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ICNTB=ICNTB+LTB(NB) 12310 

G© T® 130 12320 

770 I PAGE = I PAGE+1 12330 

WRITE (3*140) IPAGE 12340 

WRITE (3.780) IRUN 12350 

780 F0RMAT(1H .43X.32HC0MPUTED VALUES ®F X13 F0R PASS *I2//1H .5X.3HX112360 
11 *23X*3HX12*23X*6HG(X11) »2 1 X , 6HG ( XI 2 ) » 22X , 3HX 1 3 ) 12370 

12380 

FIND XI 3 12390 

12400 

L INE=4 12410 

I CNT A = 0 12420 

I CNTB=0 12430 

I CNT C = 0 12440 

NC=1 12450 

790 I A= 1 12460 

800 I AC® I A+ 1 CNT A 12470 

IB«1 12480 

810 I BC= I B+ I CNTB 12490 

XC(NC) = (XB(NB)*GA( I A ) -XA < NA ) *GB ( IB) ) / ( GA < I A ) -GB < IB) ) 12500 

I F ( ABS< l.-XB<NB)/XC(NC) )-. 000000 1 ) 840 » 820 , 82 0 12510 

820 I F ( lSL3)830*830»870 12520 

830 NC=NC+1 12530 

G® T© 870 12540 

840 I SL3= 1 12550 

I F ( ISL4 ) 850*850*860 12560 

850 I SL4= 1 12570 

XC(1)=XC(NC) 12575 

NC=2 12580 

G© T© 870 12590 

860 NC=NC+1 12600 

870 IF(NC-2000)880*740»740 12610 

880 IF( IB-LTB(NB) >890,900*900 12620 

890 IB=IB+1 12630 

G© T© 810 12640 

900 IF( IA-LTA(NA) 1910,920,920 12650 

910 I A= I A+l 12660 

G© T© 800 12670 

920 NC=NC-1 12675 

IF(LINE-50)950,930,930 12680 

930 I PAGE = I P AGE+ 1 12690 

WRITE (3,140) IPAGE 12700 

WRITE (3,940) IRUN 12710 

940 F©RMAT(1H , 37X , 32HC0MPUTED VALUES ©F X13 F0R PASS »I2»12H - C0NTIN1272O 
1UED//1H ,5X,3HX11,23X»3HX12»23X»6HG(*11 ) , 2 1 X , 6HG ( X 12 ) »22X»3HX13) 12730 

G® T© 960 12740 
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950 L I NE=LINE+2 12750 

960 WRITF (3*970) X A < N A ) * XB ( NB ) » GA < I CNT A+l ) * GB < I CNTB+1 ) , XC ( I CNTC + 1 ) 12760 

970 F©RMAT( 1H0,2 (F14.7, 12X) ,2 (E15.8*12X) »E15.8 ) 12770 

IF(LTB(NB)-1 >975*975*977 12775 

973 ND=2 12777 

I A=2 12780 

I B= 1 12783 

G© T® 980 12785 

975 IF(LTA(NA)-1 ) 1070,1070,973 12790 

977 ND=2 12793 

I A = 1 12795 

I B = 2 12800 

980 I AC= I A+ ICNTA 12810 

990 L I NF=LI NE+1 12820 

I BC= I B+ 1 CNT B 12830 

NDC=ND+ I CNTC 12840 

IF(LINE-56)l0l0»l0lC,lb00 12850 

1000 IPAGF=IPAGE+1 12860 

WRITF (3,140) I PAGE 12870 

WRITF (3*940) IRUN 12880 

L I NE=6 12890 

WRITF (3*970) X A ( NA ) * XB ( NB ) * GA ( I AC ) * GB ( I BC ) * XC < NDC ) 12900 

G® T® 1030 12910 

1010 WRITF (3* 1020) GA ( I AC ) . GB ( I BC ) * XC ( NDC ) 12920 

1020 F©RMAT(1H »52X,2(E15.8*12X) ,E15.8) 12930 

1030 IF( IB-LTB(NB) >1040,1050,1050 12940 

1040 I R= I R+l 12950 

ND=ND+1 12960 

G© T® 990 12970 

1050 I F ( IA-LTA(NA) >1060,1070*1070 12980 

1060 I A= I A+l 12990 

ND=ND+ 1 13000 

G© T© 980 13010 

1070 I CNT C = NC 13020 

1080 IF(NB-LMTB)l090, 1100, 1100 13030 

1090 ICNTB=ICNTB+LTB (NB) 13040 

NR=NB+1 13050 

G© T© 790 13060 

1100 I F ( ISL5)1110, 1110, 1120 13070 

1110 I SL 5 = 1 13080 

LMTG= ICNTB+LTB (NB) 13090 

1120 IF(NA-LMTA)1130, 1140, 1140 13100 

1130 ICNTA= ICNTA+1 13110 

N A=NA+1 13120 

NB= 1 13130 

G® T© 790 13140 
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Appendix C 
FORTRAN SYMBOLS 


Fortran Math 

Symbol Symbol Description 


A b 


Current trial solution for S(I) 


A1 

Vi 

A2 

Vi 

A60 

A 60 

A61 

V 

A62 

V 

A63 

A 63 

A64 

VO 

< 

A65 

V 

A66 

A 

66 

AD 

a 

B1 

Vi 

B2 

Vi 

CO 

c 0 

Cl 

C 1 

C2 

C 2 

CHECK 

DO 

d 0 

D1 

d , 

D2 

d 2 

D3 

DELTAX. 

d 3 


> 





Intermediate values in finding y and z 


Coefficients of F(s) 


Previous trail solution for S(I) 
Intermediate values in finding y and z 


Intermediate values used in finding 
the coefficients of F(s) 

Value of left-hand side of equations 
in cluster 1 


Intermediate values used in finding 
the coefficients of F(s) 


Increment applied to x x when finding 
suitable trial values 


Equation 

Number 

3,9 

12 


4 

3,9 

11 

5 
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Fortran 

Symbol 

Math 

Symbol 

Description 

Equation 

Number 

EO 

El 

E2 

E3 

■> 

e o 

e i 

e 2 

^ , 


Intermediate values used in finding 
the coefficients of F(s) 

7 

F3 

F5 

F7 

F9 

Fll 

f 3 , i ] 

f 5 , i 
f 7 , i 
f 9,i 

f H,i ^ 

, 

Derivatives of the Legendre poly- 
nomials, P. (x) , j = 3, 5, 7, 9, 11 

2 

FA 

F(b) 


The current value of F(s) 

3,9 

FD 

F(a) 


The previous value of F(s) 

3,9 

GA 

GB 

g( X lA>l 

g( x i B ) J 


Used in making the value of x 
satisfy the fifth equation in (1) 

16,17 

H 

F<s n )/F'(s„) 

Iterative value subs trac ted from s in 
finding the solution of F(s) 

8 

I 



Running index 


IA 



Running index for GA in generating 
new Xj 


IAC 



Actual location of GA in block 


IB 



Running index for GB in generating 
new Xj 


IBC 



Actual location in GB in block 


IC 



Loop counter in finding solution 
for F(s) (up to 10) 


ICNTA 



Lower bound for index of XA for 
given NA and NB 
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Fortran 

Math 


Equation 

Symbol 

Symbol 

Description 

Number 

LINE 


Line count on page (used in 
print control) 


LMTA 


Total number of x 1A (or XA) 


LMTB 


Total number of x JB (or XB) 


LMTC 


Total number of new x 1 (or XC) 


LMTG 


Total number of values of GB 
generated 


LT 


Number of values of S or GB for a 
given x u or x J2 (4uses) 


LTA 


Number of values of S or GB for a 
given x n (2uses) 


LTB 


Number of values of S or GB for a 
given x 12 (2 uses) 


LTD 


Dummy limit for D0 Statement (to 
eliminate subscript) 


MAXRUN 


The maximum number of passes 
allowed in finding new Xj (read in) 


NA 


Running index for XA 


NB 


Running index for XB 


NC 


Index to XC (new x t ) 


ND 


Running index for XC (new x x ) 


NDC 


Actual location of XC (new x x ) 


R2 

r 2 

/Solved solutions with satisfac-1 

18 

R3 

r 3 

ItoryXj / 

19 
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Fortran 

Symbol 


Math 

Symbol 


Description 


Equation 

Number 


ICNTB 

ICNTC 

IPAGE 

IRUN 

ISL1 

ISL2 

ISL3 


ISL4 


ISL5 


IT 

J 

J2 

J3 

JB 


Lower bound for index of XB for 
given NA and NB 

Lower bound for index of XC(NC) for 
a given NA and NB in printing 

Page number (used in printing) 

The pass number (up to MAXRUN) 

~\ 

0: G u branch 

1: G 12 branch 

2: end G branch 
12 

0: finding G t 
1: finding G 12 

2: making final computation 

0: No satisfactory new has 
been found 

1: A satisfactory new x x has 
been found 

0: Normal 

1: Counter reset when first 
satisfactory new x x has been 
found 

0: First time around in finding 
new Xj 

1: Additional times around in 
finding new x t - don't re- 
calculate LMTG 

Index for CHECK subscript 

Index (general) 

j 2 f Solved solutions with satisfac- "1 

j 3 \tory xj J 

Running index for GB 
25 


sense 

light 

simulators 


20 
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n n n n n r> 


1140 

LMTC=NC 


13150 


I F( I RUN-MAXRUN) 1150.1200*1200 

13160 

1150 

I RUN= I RUN+1 


13170 




13180 


SHIFT VALUES ©VER F©R NEXT 

PASS 

13190 




13200 


D© 1160 NB=1 »LMTB 


13210 


XA(NB)=XB(NB) 


13220 

1160 

LTA ( NB ) =LTB ( NB ) 


13230 


D© 1170 IBC=1.LMTG 


13240 

1170 

6A ( I BC ) =GB ( IBC) 


13250 


I SL5=0 


13260 


D© 1180 NC= 1 » LMTC 


13270 

1180 

XB ( NC ) =XC ( NC ) 


13280 


LMTB=LMTC 


13290 


I F ( ISL3)120. 120.119.0 


13300 

1190 

I SL2=2 


13310 


G© T© 120 


13320 




13330 


ERR0R PRINT0UTS AND ENDINGS 

13340 




13350 

1200 

WRITE (3.1210) MAXRUN 


13360 

1210 

F©RMAT ( 1H0.34HTHIS PR0BLEM 

WILL N0T C0NVERGE IN .I2.11H I TERAT I©NS 13370 

1///) 


13380 


G© T© 713 


13390 

1220 

WRITE (3.1230) 


13400 

1230 

F©RMAT ( 1H0.40HTHIS ©CURRED 

AT THE START ©F THE PR0GRAM///) 

13410 

1240 

RETURN 


13420 


END 


13430 
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Fortran 

Symbol 


Math 

Symbol 


Description 


Equation 

Number 


s 

s. 

1 

Solution of F( s) 


SDUM 


Dummy value used in computing 
CHECK 


T 

t. 

1 

Solution of equation cluster 1 

10 

TDUM 


Dummy value used in computing 
CHECK 


XI 

X 1 

Solution of equation cluster 1 (also, 
initial value of trial x x (read in) ) 


X2 

X 2 

Given solution of equation cluster 1 
(read in) 


X3 

X 3 

Given solution of equation cluster 1 
(read in) 


XA 

X 1A 

First acceptable trial x x 


XB 

X 1B 

Second acceptable trial x x 


xc 

X 1 

New value of Xj 


XDUM 

X 2 

1 

Dummy value used in computing 
F3-F11 


XINT 


Dummy x x used in computation 


XMAX 


Maximum allowable for x 1 (read in) 


Y 


Solution of equation cluster 1 

15 

YOZ 

(y/ z )i 

Intermediate value in finding Y and Z 

13 

Z 

Z. 

i 

Solution of equation cluster 1 

14 
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Appendix D 
ERROR PRINTOUTS 


As an aid to the user, a series of error printouts has been incorporated in 
the program to indicate where the problem has gone astray. A description of 
these printouts follow. 

Statement 

Number Test of Printout and its Meaning 


353 THIS PROBLEM CANNOT BE SOLVED WITH THE GIVEN 

VALUES 

No satisfactory solution can be found with the range of solutions 
that are allowed when this occurs, the program terminates. 


360 ( X j) IS NOT A SATISFACTORY VALUE 

and 604 

The current trial value of X will not work because no solution 
of F(s) can be found. The program tries another X. 

420 (s.) CONVERGES TOO SLOWLY WHEN XI 1= (x,,) AND 

X12= (x 12 ) 

The value of s in question will not converge within ten iterations 
when using the Newton-Rapson method of finding the solutions 
of F(s) . Another value of S is tried. 


590 NO POSITIVE T AVAILABLE 

None of the solutions of F(s) will give a positive value of T. 
Therefore our x t is not a satisfactory value, and another x ; is 
tried. 
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Statement 

Number 

740 


1200 


1220 


Test of Printout and its Meaning 
HELP— I AM BEING SQUEEZED OUT 


An excess of 2000 values of x n , x J2 , x J3 , g(x n )or g(x 12 ) occurs. 
When this happens, the program terminates 


THIS PROBLEM WILL NOT CONVERGE IN (MAXRUN) 
ITERATIONS 

A satisfactory convergence of x J3 will not occur for a given part 
of Xj's. Another x x is tried. 


THIS OCCURRED AT THE START OF THE PROGRAM 

The error encountered occurred before any satisfactory solution 
of x J3 was found. The program terminates. 
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Appendix E 
SAMPLE PROBLEM 


For a test run, the program was run with the following input data: 

X2 = .5917002 
X3 = .8710000 
MAXRUN = 10 

XI = .2000000 
XMAX = .2150000 
DELTAX = .001 

The first two acceptable values of XI were .2090000 and .2110000. 

When the iterative interpolation (or extrapolation) was applied to these values 
of XI, we get 


Iteration 

1 

2 

3 

4 

5 

6 
7 


with 

XI = 
X2 = 
X3 = 

we solved for R2, R3, J2, J3, and got 

R2 = 
R2 = 
J2 = 
J3 = 


New XI 

.2095635 

.2096869 

.2098787 

.2098228 

.2098287 

.2098290 

.2098290 


.2098290 

.5917002 

.8710000, 


1.003553 

.998336 

0 . 

0 . 
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Substituting these solutions into the five equations, we obtain 


Equation 1: 
Equation 2: 
Equation 3: 
Equation 4: 
Equation 5: 


-.1117587 x 10 -7 
.2607703 x 10- 7 
.3725290 x 10" 8 
.2980232 x 10“ 6 
.1546740 x 10- 3 


Since equation 5 contains fourth-degree terms, any error will be exaggerated 
when checking. 
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